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1.  Introduction 


The  primary  motivation  for  the  current  study  comes  from  our  recent  work1  on  wave 
propagation  in  microcracked  media  under  prestress;  our  interest  in  this  report  is  to 
gain  an  understanding  of  why  compressive  stresses  of  ayy  =  —90  MPa  develop 
between  stationary  (non-propagating)  fractures  in  a  medium  under  a  remote  ten¬ 
sile  prestress  of  <ryy  =  50  MPa  (cf.  Fig.  11(a)  in  Sahane  et  al.1  reproduced  in  part 
here  in  Fig.  1).  Because  of  the  stress  concentrations  at  the  crack  tips,  the  maximum 
tensile  stress  ayy  =  490  MPa  is  an  order  of  magnitude  higher  than  the  applied 
tensile  far-held  prestress,  but  closely  arranged  parallel  fractures  areas  of  compres¬ 
sive  stress  (dark  blue)  can  develop  that  are  significant  in  magnitude  relative  to  the 
applied  tensile  prestress. 


Fig.  1  An  Abaqus2  simulation1  showing  that  a  fractured  medium  under  a  tensile  far-field  pre¬ 
stress  of  50  MPa  can  generate  compressive  stresses  (dark  blue  coloration)  between  fractures 
in  close  proximity 

The  problem  of  wave  propagation  in  a  prestressed  medium  was  originally  studied 
by  Biot3'4  and  later  applied  to  problems  in  geophysics  related  to  self-gravitation 
of  the  spherical  earth,5  subsurface  detonations,6  the  propagation  of  Rayleigh7  and 
Scholte8  waves,  and  borehole  diagnostics.9  Applications  that  involve  prestressed 
media  range  from  “compliant”  materials,  (e.g.,  vascular  materials,10  dielectric  poly¬ 
mers11  and  magnetorheologic  elastomers12),  to  relatively  “brittle”  materials,  (e.g.,  con¬ 
crete,13  ceramics,14  and  glass15). 
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The  complexity  of  the  prestress  mechanics  problem  increases  dramatically  when 
one  considers  the  influence  of  microcracks  on  wave  propagation  and  other  proper¬ 
ties,  such  as  electrical  properties;  such  problems  can  be  addressed  using  the  class  of 
homogenization  methods  known  as  generalized  self-consistent  methods  (GSCMs). 16-19 
GSCMs  permit  derivation  of  the  anisotropic  effective  moduli  for  the  medium  in 
cases  where  the  cracks  are  relatively  dilute,  i.e.,  are  non-interacting  and  stationary. 
The  problem  in  which  discrete  cracks  interact  with  each  other,  but  are  stationary, 
was  addressed  by  Sahane  et  al.,1  and  the  solution  of  the  more  general  case  where 
the  microcracks  both  interact  and  propagate  across  spatial  scales  requires  the  use  of 
concurrent  multiscale  methods.20* 

To  answer  the  question  posed  earlier,  in  this  report,  we  study  the  mechanics  of  the  2 
overlapping  parallel  crack  problem  under  remote  tension,  and  show  that  this  prob¬ 
lem  is  mechanically  equivalent  to  the  problem  of  a  single  crack  parallel  to  a  rigid 
boundary  (shear-stress  free  symmetry  plane)  under  remote  tension  (see  Fig.  2(c)). 

In  Section  2,  we  present  the  method  of  integral  transforms21  to  solve  the  problem 
of  a  crack  parallel  to  a  rigid  boundary  under  remote  tension.  In  Section  3,  we  apply 
the  methods  developed  by  Erdogan  and  Gupta22-23  and  develop  a  system  of  singu¬ 
lar  integral  equations  of  the  first  kind  (with  a  Cauchy-type  singularity),  which  is 
numerically  solved  using  Gauss-Chebyshev  integration.  Fortunately,  stress  inten¬ 
sity  factors  (SIFs)  for  this  problem  are  derived  using  the  Schwartz-Neumann  alter¬ 
nating  technique24  (or  method  of  successive  approximations)  by  Chang  and  Ma25 
who  provide  a  table  of  SIFs  which  they  compare  with  the  earlier  work  of  Yoko- 
bori  et  al.26  and  Kamei  and  Yokobori.27  The  SIFs  we  derive  in  Section  4  of  this 
report  are  in  excellent  agreement  with  the  earlier  works  that  use  different  meth¬ 
ods  for  their  derivation.25-27  In  Section  5,  we  derive  expressions  for  the  normal 
stress  a yy(x,y)  and  shear  stress  crxy(x,  if)  fields  in  the  region  between  the  crack 
and  the  rigid  boundary  (symmetry  plane)  and  demonstrate  that  this  region  is  in¬ 
deed  in  a  state  of  compression,  as  our  prior  finite  element  [FE]  simulations  indicate 
in  Fig.  1,  despite  the  remotely  applied  tensile  prestress.  We  also  demonstrate  that 
both  the  SIFs  and  stress  fields  derived  via  numerical  solution  of  the  singular  integral 
equations,  compare  well  with  those  determined  using  the  commercially  available 
Abaqus/Standard2  FE  code.  Section  6  focuses  on  estimating  the  relative  error  in 

*The  commercial  multiscale  FE  code,  MultiMech,  http://multimechanics.com/  was  developed 
in  part  under  US  Army  Research  Laboratory  contract  No.  W91 1NF-07-D-  0001. 
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our  numerically  derived  solutions  by  specializing  them  to  the  closed-form  analyt¬ 
ical  solutions  found  in  Zehnder;28  our  numerical  solutions  are  found  to  be  within 
machine  precision  of  the  closed-form  analytical  solutions.  In  Section  7  we  use  this 
result  to  estimate  the  relative  error  between  the  normal  and  shear  stress  distribu¬ 
tions  using  the  numerical  solution  of  the  singular  integral  equations  and  Abaqus  FE 
results.  Section  8  outlines  the  Abaqus  FE  model  geometry  and  boundary  conditions 
used  to  solve  the  boundary  value  problem;  the  crack  is  modeled  in  Abaqus  using 
the  extended  Finite  Element  Method  (XFEM),  which  is  used  to  predict  FE-based 
SIFs,  which  are  compared  with  those  determined  using  our  numerical  solutions. 
Conclusions  follow  in  Section  9. 

2.  Governing  Equations 

This  section  outlines  the  solution  method  for  the  crack  boundary  value  problem 
illustrated  in  Fig.  2(c),  where  the  crack  interval  (—a,  +a),  without  loss  of  generality, 
can  been  normalized  to  (—1,  +1).  This  solution  can  be  obtained  by  superposition 
of  the  solution  to  the  problem  of  a  uniformly  pressurized  crack  shown  in  Fig.  3(b), 
with  the  solution  to  the  problem  of  a  semi-infinite  medium,  without  a  crack,  but 
uniformly  loaded  at  infinity  ayy  =  a_00  =  P0  shown  in  Fig.  3(c). 


Oyy  <7oo 


Fig.  2  A  fractured  medium  subjected  to  a  remote  tensile  load;  (a)  Abaqus2  simulation,1  and  the 
equivalence  of  the  (b)  2  overlapping  parallel  crack  problem  with  the  (c)  single  crack  parallel 
to  a  rigid  boundary  (shear-stress  free  symmetry  plane)  problem 


Approved  for  public  release;  distribution  is  unlimited. 


3 


(a)  (b)  (c) 


Fig.  3  (a)  The  solution  to  the  problem  depicted  in  Fig.  2(c)  is  obtained  by  superposition  of 
(b)  the  solution  to  the  problem  of  a  uniformly  pressurized  crack  Fig.  3(b),  and  (c)  the  solution 
to  the  problem  of  a  semi-infinite  medium,  without  a  crack,  but  uniformly  loaded  at  infinity 

G yy  &—oo  P0  Fig.  3(C) 


Inasmuch  as  the  solution  to  the  problem  depicted  in  Fig.  3(c)  is  trivially  uniform 
(i.e.,  ayy  =  a_oo  =  P0  everywhere,  since  the  fixed  boundary  is  shear  stress  free), 
this  section  will  focus  on  the  solution  to  the  problem  shown  in  Fig.  3(b)  by  special¬ 
izing  the  methods  and  nomenclature  outlined  in  Erdogan  and  Gupta22  for  flaws  in 
multilayered  media.  For  a  uniformly  pressurized  crack  (Jyy(x ,  0)  =  P0  on  M  <  k 
the  appropriate  boundary  conditions  on  the  crack  line  y  =  0,  are 


yy^X,  0)  — 

-Po(x)  =  Po  ; 

\x |  C  1 , 

xyiXi  0) 

0  ; 

— OO  <  X  <  oo 

v(x,  0)  = 

0  ; 

lxl  >  1 , 

and  all  components  of  stress  and  displacement  vanish  as  \J x2  +  y2  — »  —  oo  (Fig.  3(b)). 
It  is  further  assumed  that  the  medium  is  in  a  state  of  plane  strain  and  because  of  ge¬ 
ometric  and  load  symmetry  the  displacements  can  be  written  as  half-range  Fourier 
transforms  of  undetermined  functions  $  and  \D,  that  is, 


2 

ui(x,  y)  —  —  $i(£,y)sm(£x)d£, 

71  Jo 

2 

Vi(x,  y)  =  -  Vi(£,  y )  cos(£x)  d£  , 

n  Jo 


(2) 
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and  here,  i  =  —1  or  1,  where  the  subscript  i  =  —  1  refers  to  fields  in  the  infinite 
region  below  the  crack  line  y  —  0,  and  the  subscript  i  —  1  to  fields  in  the  region 
between  the  crack  line  y  —  0  and  the  rigid  boundary  y  =  h  illustrated  in  Fig.  3(b). 
On  taking  the  appropriate  partial  derivatives  of  Ui(x,y )  and  Vi(x,y )  in  Eq.  2  and 
substituting  them  into  the  2-dimensional  plane  strain,  Cauchy-Navier  equations  for 
an  isotropic  elastic  medium,  viz., 


t  \  ,o  ^d2Ui{x,y)  d  2Ui(x,y)  d\(x,y) 

(Aj  +  2 fj,i) - - - b  Hi - - h  (Aj  +  fXi) — Q  Q  —  —  0  , 


dx2 


dy 2 


dxdy 


(3) 


d\(x,y )  ,  0  d\(x,y)  ,  ^d2Ui(x,y)  _ 

^ i — - b  (Aj  +  2 in) - 7—7 - b  (A i  +  jdi) — 0  0  —  —  0  , 


dx2 


dy 2 


dxdy 


and  employing  Fourier’s  sine  formula, 


f(x)  =  Fs{Fs(f (x)  :  x  xj  0  <  x  <  00  ,  (4) 

(e.g.,  Eq.  1.1.5  of  Titchmarsh29)  where  Ts  denotes  the  half-range  sine  transform, 

Ts{i)  =  T s(f(x)  :  x  ->•  f)  =  y?  I  f(x)  sin(£z)  dx  ,  (5) 

(e.g.,  Eq.  1.2.3  of  Titchmarsh29)  we  arrive  at  the  following  2  ordinary  (coupled) 
differential  equations  in  <f>  and  T : 

-Ktiem,v)  +  Ka&l&y)  -  y)  =  0, 

(6) 

~Ki2e%^,y)  +  KaW&y)  +  £*<(£,3/)  =  0. 

Here  Kt]  =  =  2  —  2z/j,  and  Ka  =  =  1  —  2z/j,  where  Aj  and  /j,:  are  the 

Lame  and  shear  moduli,  is  Poisson’s  ratio,  and  the  primes, denote  derivatives 
with  respect  to  y.  On  solving  Eq.  6  for  <f>j(£,  y)  we  arrive  at  the  following  fourth- 
order  ordinary  differential  equation: 


with  general  solution, 


y)  —  (An  +  Ai2)e  ^ y  +  (Aft  +  Aii)e^y  .  (8) 
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On  substituting  Eq.  8,  and  its  second  derivative  into  Eq.  61,  and  integrating  y) 
over  all  y  results  in  the  general  solution  for  \Eq(£,  y). 


^ifay)  =  (Ai  +  (y  +  y)Ai2)e  iy  +  (-A3  +  Cj  -y)Ai4)e^y  ,  (9) 

where  n,  =  3  —  4//,.  In  these  equations,  the  Arj  =  Atj  (£) ,  i  =  —1  or  1,  andj  = 
1, ...,  4  are  to  be  determined  from  the  boundary  and  continuity  conditions  for  the 
specific  problem.  On  substituting  Eq.  8  and  Eq.  9  into  Eq.  2  we  arrive  at  the  general 
expressions  for  the  displacements,  given  without  derivation  in  Erdogan  and  Gupta,22 

2 

Ui(x,y)  =  -  {(Ai  1  +  Ai 2)e~iy  +  ( Ai3  +  AiA)e^y}  sin(£z)  d£, 

71  Jo 

(10) 

2  r°°  k  ■  k 

Vi(x,  y)  =  -  {(An  +  (-J  +  y)Ai2)e~iy  -  (Aa  -  (-J  -  y)AiA)eiv}  cos(£z)  d£  . 
77  Jo  £  4 

The  normal  alyy(x,y)  and  shear  crlxy(x,y)  stresses  follow  immediately  by  taking 
the  appropriate  derivatives  of  Eq.  10  and  substitution  into  Hooke’s  law  written  in 
Cartesian  coordinates, 


vL(x,y)  =  (Aj  +  2  Hi) 


dv.i(x,y)  dui(x,y) 


dy 


+  A 


dx 


V'xyfav) 


l^i 


d  Uj(x,y)  d  Vj(x,y) 

dy  dx 


(ID 


which  results  in 


a. 


yy 


A,y) 


2/ij 


2 

71 


[£(Alii  +  Ai2y)  +  2(1  —  Vj)Ai2\e  ^ y 


+[-£(A3  +  Auy)  +  2(1  -  Ui)AiA]eiy}  cos(£x)  dg, 

AtlA  =  1  +  Aay)  +  (1  -  2 Vi)Aa]e^ 

2/ij  77  J 0 


(12) 


+  [£(A3  +  Any)  -  (1  -  2 u^A^e^}  sin(£z)  d£  . 
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3.  Boundary  and  Continuity  Conditions 


The  geometry  for  this  problem  can  be  envisioned  as  consisting  of  2  simply  con¬ 
nected  regions:  one  region  consists  of  an  isotropic  elastic  strip  of  width  h,  where 
0  <  y  <  h  and  —  oo  <  x  <  oo,  that  is  bonded  to  another  region  that  consists  of  a 
semi-infinite  isotropic  elastic  plane,  where  —  oo  <  y  <  0  and  — oo  <  x  <  oo.  The 
bond  connecting  the  first  region  to  the  second  region  occurs  along  y  —  0,  and  all  x, 
except  for  \x\  <  1  (see  Fig.  3(b)).  Eight  constants  are  required  to  solve  this  prob¬ 
lem,  4  for  the  strip  and  4  for  the  semi-infinite  plane.  Two  constants,  A_u  and  A_\2 
are  identically  zero  for  the  displacements  to  be  bounded  as  y  — *  —  oo.  This  leaves  6 
constants  to  be  determined  from  the  following  6  conditions: 


1)  a1xy(x,h)  =  0  ;  y  =  h, 

2)  vi(x,  h)  =  0  ;  y  =  h, 

3)  aly(x,0)  -  a-i(x,0)  =  0  ;  y  =  0, 

(13) 

4)  crly(x,  0)  —  -cr~y(x,  0)  =  0  ;  y  =  0, 

5)  .°(ui^«-l)  =  f15(x  -  t)  ;  y  =  0, 

6)  d{viexv-l]  -  h  5(x  -  t)  ;  y  =  0. 

Equations  13i  and  132  are  the  stress  boundary  conditions  on  the  upper  strip  rigid 
boundary  depicted  in  Fig.  3(b).  Equations  133  and  134  are  the  stress  continuity 
conditions  on  the  crack  line.  Equations  13s  and  136  represent  the  symmetrically 
disposed  unit  dislocations  at  y  =  0,  x  =  t.  In  the  present  problem,  therefore,  a  6  x  6 
system  is  inverted,  which  provides  expressions  for  An,  A12,  A13,  Au,  A_13,  A_14 
as  functions  of  fi,  f2,  and  The  system  of  equations  are  derived  by  substitution 
of  Eq.  10  and  Eq.  12  into  Eq.  13,  the  result  of  which  appears  in  matrix-vector  form 
as 

Ax  =  b  ,  (14) 
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where 


A  = 


e-** 

e“^  (h  +  |j 

—eh^ 

eK 

0 

0 

■e-^e 

e~h^(2u  —  —  1) 

eh^ 

eh^(2v  +  h£  —  1) 

0 

0 

-2(1  -  i/) 

-i 

2(1  -  u) 

£ 

-2(1  -  v 

2/y  -  1 

i 

2v  -  1 

l-2u 

0 

i 

0 

0 

—  K 

i 

—  K 

K 

,(15) 


and  we  have  dropped  the  subscript  i  in  this  equation  since  the  elastic  material  above 
and  below  the  crack  line  is  assumed,  for  this  problem,  to  be  identical  (i.e.,  v  = 

ui  =  z/_!,  and  k  =  3  —  4  v),  while 


'  Tin  ^ 

r  0 

-4l2 

0 

A\s 

0 

,  b  = 

v4i4 

0 

-4-13 

fi  cos  (ft) 

A- 14 

_  /2sin(ft)  _ 

(16) 


The  expressions  for  the  x  =  Aij( £)  are  obtained  by  finding  the  inverse  of  ma¬ 
trix  Eq.  15,  x  =  A~l  b,  viz., 


x  = 


cos(tg)/i  .  (l-2v)  sin (fg)/2 

2£  _r  2(-l+k+2v)£ 

cos(tg)/i  .  sh\(tj)f2 
-4+4 v  '  2-2k-Au 

e~2h((_l +i/+fc£)  cos(t£)/i  |  e-2hZ(l-2v-2ht)sin(tt)f2 

2(— 1 +v)£  +  2(-l+fc+2i,)£ 

e~2^  cos(tg)/!  e~2^  sin (tg)/2 

2(2-2i/)  “l"  2(-l+k+2u) 

e  ( —  l  +  ['+/l£)'\  f  1  \  c 

J  cos(t£)/i  ^  e— 2h.{ (i_|_e2/i£(i_2i/)_ 2i/— 2fe^)  sin(t£)/2 
2|  I  2(-l+fc+2i/)£ 

(-l+e_2h?)  cos(t£)/i  (l+e-2^)  sin(t£)/2 
4-4i/  '  2(-l+fc+2i/) 


•  (17) 
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A  compact  expression  for  the  representation  of  the  in  Eq.  17  is 


AAO  =  F±  (0  fi  cos (ft)  +  i^(0  h  sin (&) , 
(*  =  “!>  !»  7  = 


(18) 


where  the  F-j  are  functions  of  the  elastic  properties  and  f.  On  substituting  Eq.  18 
into  Eq.  12  gives  the  2  stress  components  ayy  and  axy  due  to  the  2  unit  dislocations 
/i  and  f2  given  by  Equations  135  and  136  in  the  upper  strip  as 


ayy(x,y,t)  =  hn(x,y,t)  fa  +  h12(x,  y,t)  f2  , 
crxy(x,  y,  t )  =  h21(x,  y,  t )  fx  +  h22(x,  y,  t )  f2  . 


(19) 


Equation  19  represents  the  Green’s  functions  for  the  current  problem  and  the  hl;} 
are  singular  integrals  involving  the  F^.  Assuming  now  that  the  unit  dislocations,  f\ 
and  f2,  are  functions  of  t  and  represent  the  unknown  distributions  of  dislocations 
on  the  crack,  the  integral  expressions  of  Eq.  19  for  the  determination  of  fi(t)  and 
f2(t )  become 

lim  /  hn(x,y,t)  fi(t)  +  h12(x,y,t)  f2(t)dt  =  ayy(x ), 
y^°  Jo 

(20) 

lim  /  h21(x,y,t)  /i(t)  +  h22(x,y,t)  f2(t)dt  =  axy(x) , 
y^°  Jo 

o  <  x  <  l , 


or 


4(i, -1) 


2(I  °yv(x,°)  I  h(t)dt  J  -2ft£cos((S  -  x)()e  2h!  d( 


1  pi  poo 

-  /  f2(t)dt  /  — (1  +  2 /if)  sin((t  —  x)£)e~2h^  df 

71  J- i  Jo 


(21) 


—  lim  —  /  f2(t)dt  /  sin((t  —  x)f)e  df  , 

v-+o  vr  J-i  Jo 
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and 


A(u-1)  .  .  1 

— — — <Jxy(x,0)  =  -- 

2j  f  t  7T 


f2(t)dt  /  — 2/r£cos((f  —  a;)£)e  2h^dt; 


'-i 


1 

7T 


I  fi(t)dt  J  (1  -  2/i£)sin((t  -x)£)e 


(22) 


1  f1 

+  lim  —  /  fi(t)dt 

y^O  7T  y_1 


-x)Oe 


Integration  of  the  sine  and  cosine  transforms  in  the  iterated  integrals  in  Eq.  21 
and  Eq.  22  results  in 

-  [  P^-dt  +  -  [  ku(x,t)  fi(t)  +  k12(x,t)  f2(t)dt  =  0, 
vr  J_xt-x  tt  J 


7 r  ,/_!  t  —  X  71 


k2i(x,t)  fi(t)  +  k22(x,t)  f2(t)dt  = 


'-i 


(k  +  1)P0 

2/i 


,  (23) 


-  1  <  x  <  1 . 

In  the  derivation  of  Eq.  23  we  have  assumed  that  the  applied  normal  stress  ayy  (x,0) 
on  the  crack  is  constant  and  the  applied  shear  stress  <rxy(x,  0)  on  the  crack  is  zero. 
In  addition,  the  unknown  dislocation  densities  /i,  /2  in  Eq.  23  are  written  as  the 
tangential  derivatives  of  the  relative  crack  displacements,  as 

,  /  \  _  d[ui(x)  -  u- i(x)]  ,  _  d[vi(x)  -  n-i(x)] 


(24) 


j  i  w 


dx 


dx 


<7xy(x,  0)  =  0  ,  <Tyy(x,  0)  =  -P0 


(!  +  «) 

2/i  ’ 


-1  <  x  <  1 . 
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The  kernel  functions  in  Eq.  23  are  symmetric  k,3  =  k:]l  and  are  derived  herein  for 
the  specific  problem  illustrated  in  Fig.  3(b)  as 


ku(x,t) 


(4 h2  —  (t  —  x)2)  (t  —  x) 
(4 h2  +  (t-  x )2f  ’ 


ku(x,t) 


k2i(x,t) 


2 h  (Ah2  —  (t  —  x )2) 
(4 h2  +  (t-x)2)2  ’ 


(25) 


k22(x,t) 


(12 h2  +  (t  —  x)2)  (t  —  x) 
(Ah2  +  (t-  x)2)2 


The  unknown  dislocation  densities  f\ (t)  and  f2(t )  are  numerically  determined  by 
reducing  the  singular  integrals  in  Eq.  23  to  an  infinite  system  of  linear  algebraic 
equations  using  the  orthogonality  properties  of  the  Chebyshev  polynomials.23  It  is 
common  practice  to  write  the  dislocation  densities  in  terms  of  an  infinite  series, 
truncated  at  the  k  =  Nth  term,  of  Chebyshev  polynomials  of  the  first  kind,* 

1  OO 

hit)  =  -7==J2AkT2k{t), 

v71  - 12  t~i 

(26) 

1  X 

f2(t)  =  ,  y>fcr2fc-i(t)  • 

V71  - 12  k=1 

On  substituting  Eq.  26  into  Eq.  23,  and  using  the  orthogonality  conditions  of  the 
Chebyshev  polynomials,  one  can  derive  a  set  of  functional  equations  involving  the 
unknown  constants  Ak  and  Bk  (cf.  Eq.  (7.97)  of  Erdogan  et  al.23).  The  functional 
equations  can  be  solved  using  a  weighted  residual  method  to  arrive  at  the  following 

*The  more  general  Jacobi  polynomials  are  used  in  solutions  to  problems  where  f(t)  = 

g(t)(l—t)a(t+l)a  in  Eq.  23  and  g(t)  =  Ykkk=o  since  the  solution  of  Eq.  23  has  integrable 

singularities  at  the  end  points,  then  «  =  /?  =  — 1/2,  and  the  Jacobi  polynomials  coincide  with  the 
Chebyshev  polynomials  of  the  first  kind  used  in  Eq.  26  (cf.  discussion  on  pages  380-38 1  of  Erdogan 
et  al.23). 
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system  equations: 


N 


n  =  1 


N 


(27) 


n  =  1 


Fik  =  0,  (k  —  1,  ...N) , 


where  the  definitions  of  constants  a^n,  bkn,  Ckn,  and  dkn  can  be  found  in  Appendix  A. 

4.  Stress  Intensity  Factors 

In  this  section,  SIFs  are  computed  and  compared  with  the  prior  results  of  Chang  and 
Ma25  who  tabulate  their  results  that  are  based  on  the  Schwarz-Neumann  alternat¬ 
ing  technique,24  along  with  those  of  Yokobori  et  al.,26  and  Kamei  and  Yokobori27 
(Table  1).  Using  the  formalism  developed  in  this  report,  the  mode  I,  k\ ,  and  mode 
II,  k2,  SIFs  can  be  written  in  a  straightforward  manner,  cf.  page  394  of  Erdogan  et 


al.,23  as 


and 


where  the  infinite  series  is  normally  truncated  to  a  finite  number  of  terms,  N,  say 
k  =  N  =  9,  and  the  Ak  and  74  are  given  in  Eq.  A-7.  One  immediately  observes 
from  Table  1  that  the  SIFs  obtained  by  the  prior  authors  using  different  solution 
methodologies  agree  quite  well  with  our  own  results*  for  the  case  N  =  9,  cf.  Ap¬ 
pendix  B  that  illustrates  SIF  convergence  for  N  =  3,  6,  9  and  various  values  of 
a/h.  In  general,  as  a/ h  0,  that  is,  as  the  crack  distance  from  the  rigid  boundary 
increases,  the  mode  II  stress  intensity  factor  approaches  zero,  k2  — >  0,  and  the  mode 

*  We  use  the  Erdogan  et  al.23  convention  for  the  shear  stress  so  that  our  Fjj  values  are  opposite  in 
sign  to  that  used  by  the  other  authors,  but  consistent  with  those  computed  using  Abaqus  in  section  8. 
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I  stress  intensity  factor  approaches  unity,  k\  — »  1*.  Since  the  SIFs  derived  by  ear¬ 
lier  authors  for  the  2  overlapping  parallel  crack  problem  under  remote  tension  are 
nearly  the  same  as  our  problem  of  a  single  crack  parallel  to  a  rigid  boundary  under 
remote  tension  Table  1,  these  2  problems  are,  therefore,  mechanically  equivalent. 

The  equation  used  to  derive  the  initial  angle  of  crack  growth  (see  Table  1)  subjected 
to  mixed-mode  loading  is  written  as  a  function  of  k\  and  /c2,  and  given  by  Eq.  1  on 
page  236  of  Erdogan31  as 

k\  sin  90  +  /c2( 3  cos  6*0  —  1)  =  0  .  (29) 

For  cracks  that  are  distant  from  the  rigid  boundary,  the  initial  angle  of  crack  growth 
(cleavage  angle)  approaches  zero,  60  — >■  0,  the  crack  propagation  is  self-similar; 
cracks  that  approach  the  rigid  boundary  propagate  at  ever  increasing  angles  (see 
Fig.  4)  relative  to  the  crack  line,  and  away  from  the  rigid  boundary  (see  Fig.  3(b)). 


Table  1  Normalized  SIFs  for  2  parallel  cracks  in  an  infinite  plane  under  remote  tension, 
where  I<i  =  Fjcry/na,  Kn  =  Fucr^/Fa,  and  k\  =  Fj;  k->  =  Fjj  are  determined  using  Eq.  28 


a/h 

jp  26,27 

Fi25 

Ff 

jp  25 

r  11 

F,id 

O 

cr 

0.01 

0.99996 

-1.8747x  10~7 

0 

0.1 

0.9963 

-0.00018 

0.02 

0.2 

0.985526 

0.9858 

0.9857 

0.0014 

-0.0014 

0.16 

0.4 

0.950826 

0.9505 

0.9505 

0.0094 

-0.0094 

1.13 

0.6 

0.908926 

0.9092 

0.9086 

0.0246 

-0.0246 

3.09 

0.8 

0.872726 

0.8722 

0.8706 

0.0431 

-0.0429 

5.61 

1.0 

0.83 1927 

0.8431 

0.8403 

0.0611 

-0.0607 

8.18 

1.25 

0.803727 

0.8166 

0.8131 

0.0803 

-0.0793 

10.94 

1.66 

0.7860 

-0.1011 

14.21 

2.0 

0.756927 

0.7734 

0.7737 

0.1166 

-0.1128 

15.95 

5.0 

0.7445 

-0.1567 

22.02 

10.0 

0.7385 

-0.1792 

24.74 

a  Fj  and  Fji  values  are  from  Table  B-l,  N  =  9. 
b  do  values  are  derived  using  Fj  and  Fjj  values  in  this  Table  1. 


*Some  authors  use  the  convention  Kj  — >  dy/Fa  as  a/h  —>  oo  but  we  use  the  normalized 
constant  k\  —>  <J\fa  as  a/h  ^  oo,  cf.  footnote  2  of  Erdogan.30 
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Fig.  4  Fi,  Fu,  and  cleavage  angles  90  vs.  a/h  from  Table  1  for  a  crack  parallel  to  a  rigid 
boundary  under  symmetric  loading 
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5.  Normal  and  Shear  Stresses 


In  this  section,  we  determine  the  normal  cryy(x,  hi)  and  shear  stress  axy(x.  h\)  at 
y  =  hi,  where  0  <  hi  <  h  and  compare  our  results  with  those  obtained  using  the 
Abaqus  commercial  FE  code  where  the  crack  is  modeled  using  the  XFEM. 

On  substituting  the  Aiy(£)  from  Eq.  17  into  Eq.  12  we  arrive  at  the  iterated  integrals 
for  the  normal  (Tyy(x,  hi)  and  shear  stress  axy(x,  hi): 


- - LVyy(x,hi) 

2/i 


,-(2h+h1)i 


,-(2h+h1)i 


(^—e2h^hi  +  e2hli(-2 h  +  hi))  ^ cos ((t  -  x)£)  d£ 


(-e2,l5(l  +  hiO  -  e2hl«(l  +  2 h£  -  hi£))  sin (£(i  -  x))  d£ 


(30) 


f2(t)dt  f  -e~(2h+hi)£  fe2hZhl  +e2h^(-2h  +  h1)]  £cos((t-x)£)dt. 
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On  integrating  the  sine  and  cosine  transforms  in  Eq.  30  we  arrive  at  the  following 
equations  for  the  normal  stress  oyy(x,  hi): 


-(«  +  !) 
2(i 


@ yy  (*£>  h\  ) 


1 
7 r 


mi 


hi  (- h\  +  (t  -  x)'2) 
(hf  +  (t  -  x)2)2 


hi  ((—2 h  +  hi)2  -  (t  -  x)2)  2 h  (-(-2 h  +  hi)2  +  (t  -  x)2)  ^ 

((—2 h  +  hi)2  +  (t  -  x)2)2  ((-2 h  +  hi)2  +  (t  —  x)2)2 

(31) 

1  f1  f  Mr  2/li  1 

J  —i  h[t){(h2  +  (t-x)2)2  +  h2  +  (t-x)2 


2(-2h  +  hif 


+ 


((-2 h  +  hi)2  +  (t  -  x)2)2  '  (-2 h  +  hi)2  +  (t-x)2 

and  shear  stress  axy(x ,  h  i )  in  terms  of  polynomials  in  x  and  t. 


}(t  —  x)dt , 


-(«  +  !) 
2(i 


axy(x,hi)  =  -  fi(t){ 


-2  h\ 


+ 


7T 


'-1 


(h\  +  (t-  x)2)2  h2  +  (t-x) 


2(-2h  +  hif 


((-2 h  +  hi)2  +  (t-  x)2)2  (-2 h  +  hi)2  +  (t-  x)2 


}(t  —  x)dt 


-  f1  f  (fU  h!  ((-2h  +  hi)2  -  (t-x)2) 
"vr  ((-2 h  +  ln)2  +  (t-x)2)2 


(32) 


^2  h  (-(-2  h  +  ^i)2  +  (t  —  x)2)  +  hi(hi  +  t  —  x)(hi  -  t  +  x)  ^ 

+  ((-2h  +  hi)2  +  (t  -  x)2)2  (/i2  +  (f  -  x)2)2 

On  substituting  the  dislocation  densities  given  by  Eq.  26  into  Eq.  31  and  Eq.  32 
together  with  the  Ak,  and  Bk  determined  by  the  methods  described  in  Appendix  A, 
and  employing  Gauss-Chebyshev  numerical  integration  using  Eq.  A-3,  we  arrive  at 
the  normal*  cryy(x,  hi)  and  shear  stress  crxy(x,  hi)  distributions  for  0  <  h  \  <  h 
along  line  0  <  x  <  5,  where  they  are  compared  with  the  Abaqus  commercial 

*Recall  that  the  normal  stress  distribution  depicted  in  Fig.  3(a)  consists  of  superposing  Eq.  31 
with  a  unit  remote  stress  at  infinity  ayy  =  <r_ oo  =  1. 
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FE  solutions  (Figs.  5  and  6).  Also,  the  normal  stress  approaches  the  far- field  unit 
applied  tensile  stress  in  Fig.  5  and  far-held  null  shear  stress  in  Fig.  6. 

In  Section  6,  we  verify  the  accuracy  of  the  Gauss-Chebyshev  solutions  depicted 
in  Figs.  5  and  6  by  showing  they  are  within  machine  precision  of  analytical  solutions 
for  the  problem  of  a  crack  in  an  infinite  plane  subjected  to  remote  tension.  This 
result  allows  us  to  accurately  estimate  the  error  between  our  solutions  and  those 
obtained  using  Abaqus.  Finally,  we  show  in  Fig.  7  that  for  all  parallel  cracks  that 
are  in  close  proximity  to  each  other,  (i.e.,  at  normalized  distances  a/h  >  0.95),  the 
stress  between  the  parallel  cracks  is  compressive  (stress  is  negative  in  compression), 
and  the  compressive  stress  increases  linearly  in  magnitude  with  normalized  distance 
a/h  (Fig.  7). 


Fig.  5  Normal  stress  ayy(x,  hi)  solution  for  the  boundary  value  problem  depicted  in  Fig.  3(a): 
Abaqus  simulation  compared  with  Eq.  31  for  h  =  0.6  and  hi  =  0.2;  with  Ai,  ...An,  Bi,  ...,  BN 
constants,  k  =  N  =  6  in  Eq.  26 
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Abaqus  —  Eq.  32 


Oxy 


Fig.  6  Shear  stress  axy(x,  hi)  solution  for  the  boundary  value  problem  depicted  in  Fig.  3(a): 
Abaqus  simulation  compared  with  Eq.  32  for  h  =  0.6  and  hi  =  0.2;  with  Ai,  ...An,  B\, ...,  BN 
constants,  k  =  N  =  6  in  Eq.  26 


—  Maximum  stress 

0"yy(x,h) 


a/h 


Fig.  7  Maximum  stress  ayy(x,h)  along  rigid  boundary  line  y  =  h  vs.  a/h;  far-field  tensile 
stress  is  ayy( x,  —  oo)  =  1  with  tension  positive,  compression  negative.  At  a/h  «  0.95  the  stress 
changes  from  tension  to  compression  and  then  increases  linearly  with  a/h. 
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6.  Verification  of  Our  Derived  Solutions 


In  this  section,  we  verify  the  accuracy  of  our  derived  solutions  and  the  Gauss- 
Chebyshev  numerical  integration  methodology  against  a  closed-form  solution  for 
the  normal  and  shear  stress  distributions  in  an  infinite  plate  containing  a  crack  sub¬ 
jected  to  a  remote  tensile  stress  ayy(x ,  oc)  =  =  1.  The  normalized  analytical 

solution  for  the  stress  fields  is  given  in  Zehnder,28 


(33) 


and 


(34) 


where  z  =  x  +  iy,  and  3ft  denotes  the  real  part,  and  3  the  imaginary  part  of  the 
quantity  in  parentheses.  These  solutions  are  plotted  in  Figs.  8  and  9.  The  accuracy  of 
our  derived  normal  stress,  Eq.  31,  and  shear  stress,  Eq.  32,  solutions  is  determined 
by  specializing  the  solutions  to  the  case  when  the  crack  is  distant  from  the  rigid 
boundary  (e.g.,  at  y  —  h  =  1000),  but  the  stress  distribution  is  evaluated  close  to 
the  crack  line  y  —  0,  (i.e.,  at  a  distance  y  =  hi  =  0.1).  We  compute  the  relative 
error32  between  our  Gauss-Chebyshev  numerical  solution  and  the  exact  analytical 
solution28  for  the  case  of  a  crack  in  an  infinite  plane  subjected  to  a  remote  tensile 
stress.  If  E(/c,  hi)  represents  the  exact  value  of  the  stress  given  by  Eq.  33  or  Eq.  34 
and  T,(x,hi)  represents  the  approximate  value  of  the  corresponding  stress  using 
Gauss-Chebyshev  numerical  integration,  then  the  relative  error  is  given  as 


relerr  =  1  —  E(x,  hi)/T,(x,  hi) . 


(35) 


Illustrations  of  Log10(|relerr|),  where  |relerr|  denotes  the  absolute  value  of  relerr,  for 
the  normal  and  shear  stress  distributions  appear  in  Figs.  10  and  11  with  oscillatory 
solution  accuracy  ranging  from  10  to  20  significant  digits  using  n  =  100  terms 
and  averaging  about  15  significant  digits  using  n  =  300  terms  in  the  numerical 
integration  Eq.  A-3.  The  relative  error  is  defined  only  if  T,(x,hi)  ^  0,  since  if 
T,(x,hi)  =  0  the  error  becomes  unbounded;  a  loss  in  precision  is  indicated  by 
the  arrow  in  Fig.  11  (for  n  =  300)  at  an  ./'-position  where  the  shear  stress  passes 
through  a  zero  point  on  the  ordinate  in  Fig.  9.  The  relative  error  for  the  shear  stress 
is  rescaled  and  plotted  for  the  n  =  300  case  in  Fig.  12.  It  is  clear  that  the  relative 
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error  between  our  derived  normal  stress,  Eq.  31,  and  shear  stress,  Eq.  32,  solutions 
and  the  exact  analytical  solutions  is  at  machine  precision  using  double  precision 
arithmetic  in  Mathematica33  with  approximately  15  significant  digits  in  the  solution, 
cf.  pg.  253,  Table  7.1  of  Oberkampf  and  Roy.34 


Eq.  33 


Fig.  8  Exact  normal  stress  crvy(x,  y )  distribution  at  y  =  0.1  (Eq.  33)  in  an  infinite  plate  con¬ 
taining  a  crack  subjected  to  a  remote  stress  at  y  =  ±oo 


Approved  for  public  release;  distribution  is  unlimited. 


20 


Eq.  34 


crx y 


Fig.  9  Exact  shear  stress  axy(x,  y)  distribution  at  y  =  0.1  (Eq.  34)  in  an  infinite  plate  contain¬ 
ing  a  crack  subjected  to  a  remote  stress  at  y  =  ±oo 


Log10(|relerr|) 


n  =  100 


n  =  300 


Fig.  10  Numerical  accuracy  of  the  normal  stress  determined  with  Gauss-Chebyshev  numerical 
integration  relative  to  the  exact  solution  shown  in  Fig.  8  over  the  space  interval  0  <  x  < 
2;  over  this  interval,  the  numerical  solution  accuracy  ranges  from  12  to  18  significant  digits 
for  n  =  100  terms  and  averages  15  significant  digits  for  n  =  300  terms  in  Eq.  A-3;  with 

Ai,  ...Am,  Bi,  Bn  constants,  k  =  N  =  9  in  Eq.  26 
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n  =  100 


n  =  300 


Log10(|relerr|) 


Fig.  11  Numerical  accuracy  of  the  shear  stress  determined  with  Gauss-Chebyshev  numerical 
integration  relative  to  the  exact  solution  shown  in  Fig.  9  over  the  space  interval  0  <  x  < 
2;  over  this  interval,  the  numerical  solution  accuracy  ranges  from  10  to  20  significant  digits 
for  n  =  100  terms  and  averages  15  significant  digits  for  n  =  300  terms  in  Eq.  A-3;  with 

Ai,  ...An,  Bi,  ...,  Bn  constants,  k  =  N  =  9  in  Eq.  26 


Log10(|relerr|) 
—  13.8  r 

- 14.0 

- 14.2 

-14.4 

- 14.6 

- 14.8  'L - 

— 15.0 


n  =  300 


0.5  1.0 


Fig.  12  Rescaled  numerical  accuracy  of  the  shear  stress  for  n  =  300  terms  shown  in  Fig.  11 
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7.  Relative  Error  in  the  Abaqus  Simulations 


Earlier  in  Section  6,  we  demonstrated  that  our  numerical  solutions  are  at  machine 
precision  using  double  precision  arithmetic  in  Mathematical33  In  this  section  we 
will  evaluate  the  accuracy  of  the  normal  and  shear  stress  simulations  illustrated 
in  Figs.  5  and  6  that  were  determined  using  the  Abaqus  FE  code.  This  can  ac¬ 
complished  by  evaluating  Eq.  31  and  Eq.  32  at  the  same  (x,y  =  0.2)  coordinate 
locations  as  those  used  in  the  Abaqus  simulations,  and  calculating  the  relative  error 
using  Eq.  35. 

Figure  13  shows  the  relative  error  in  the  Abaqus  simulation  over  the  interval  0  < 
x  <  5  where  the  numerical  solution  accuracy  approaches  6  significant  digits  for  n  = 
300  terms  in  Eq.  A-3;  with  Ai,  ...An,  B\, ...,  Bn  constants,  k  =  N  =  6  in  Eq.  26. 
The  relative  error  is  greatest  over  the  interval  0  <  x  <  0.5  and  becomes  undefined 
near  x  =  0.6  where  the  normal  stress  solution  ara(0.6,  0.2)  — >-  0  (Eq.  31),  and 
therefore  the  relative  error  at  this  location  is  undefined.  Figure  14  shows  the  relative 
error  in  the  Abaqus  simulation  over  the  interval  0  <  x  <  5  where  the  numerical 
solution  accuracy  approaches  6  significant  digits  for  n  =  300  terms  in  Eq.  A-3; 
with  A\,  ...An,  B\, ...,  Bn  constants,  k  =  N  =  6  in  Eq.  26.  The  relative  error 
is  at  a  minimum  over  the  interval  0  <  x  <  2  and  becomes  undefined  near  x  = 
2.65  where  the  shear  stress  solution  axy(2.65,  0.2)  — >■  0  (Eq.  32),  and  therefore  the 
relative  error  at  this  location  is  undefined.  Figures  13  and  14  demonstrate  that  the 
XFEM  implementation  in  Abaqus  can  be  used  to  accurately  predict  the  classical 
static  stress  fields  in  a  linear  elastic  medium  with  a  crack. 
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Fig.  13  Relative  error  of  the  normal  stress  determined  with  Abaqus  (refined  mesh  model  of 
Section  8)  vs.  Gauss-Chebyshev  numerical  integration  shown  in  Fig.  5  over  the  space  interval 
0  <  x  <  5.  The  arrow  indicates  where  then  normal  stress  has  a  zero-crossing  E(.i\  h  t )  = 
(7yy(x.  0.2)  — ►  0  in  Fig.  5  and  increases  the  error  according  to  Eq.  35. 


Fig.  14  Relative  error  of  the  shear  stress  determined  with  Abaqus  (refined  mesh  model  of 
Section  8)  vs.  Gauss-Chebyshev  numerical  integration  shown  in  Fig.  6  over  the  space  inter¬ 
val  0  <  x  <  5.  The  arrows  indicate  where  the  shear  stress  has  a  zero-crossing  E(x,  hi)  = 
axy(x,  0.2)  — >  0  in  Fig.  6  and  increases  the  error  according  to  Eq.  35. 
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8.  Abaqus  Finite  Element  Modeling 


The  computational  FE  modeling  is  done  with  Abaqus/Standard  using  the  XFEM 
capability.  The  XFEM  enriches  elements  with  added  degrees  of  freedom  that  al¬ 
low  for  the  modeling  of  discontinuities,  specifically  the  discontinuous  displacement 
field  across  the  crack  surfaces  and  the  asymptotic  crack-tip  stress  field.  SIFs  can  be 
calculated  from  contour  integration  (J-integrals)  around  the  crack  tip.  While  the 
XFEM  does  not  require  the  crack  surface  or  tip  be  congruent  with  the  mesh,  be¬ 
cause  of  how  the  contour  paths  for  the  J-integral  are  computed,  a  regular  mesh  that 
is  relatively  aligned  (e.g.,  parallel)  to  the  crack  reduces  error  and  is  desirable.35,36 

The  numerical  solution  to  the  singular  integral  equations  assumes  plane  strain  con¬ 
ditions,  but  the  XFEM,  which  incorporates  asymptotic  fields  in  Abaqus,  is  only  im¬ 
plemented  in  3-dimensional  (3-D)  analyses  for  stationary  cracks.  Therefore,  bound¬ 
ary  conditions  have  to  be  applied  to  the  3-D  FE  geometry  to  simulate  plane  strain 
conditions.  The  geometry  and  boundary  conditions  are  shown  in  Fig.  15.  The  ge¬ 
ometry  uses  half-symmetry  at  x  =  0  with  the  ^-displacement  fixed  and  the  crack 
parallel  to  the  tc-axis.  The  //-displacements  are  fixed  on  the  top  surface  and  a  tensile 
distributed  load  is  applied  to  the  bottom  surface,  similar  to  Fig.  3(c).  To  approxi¬ 
mate  plane  strain  conditions,  the  ^-direction  thickness  is  small  compared  the  lateral 
dimensions  (L  =  5,  W  =  10,  T  =  1)  and  the  displacements  in  the  ^-direction 
are  fixed  on  both  ^-direction  faces  to  restrict  the  thickness  from  changing  because 
of  the  Poisson  effect.  The  mesh  is  biased  in  the  x-  and  //-directions  so  that  it  is  uni¬ 
form  and  refined  at  the  crack  tip,  with  the  bias  allowing  larger  elements  at  the  edges 
of  the  model  furthest  from  the  crack  tip.  It  will  be  shown  that  while  acceptable  re¬ 
sults  for  the  SIFs  can  be  calculated  using  a  relatively  coarse  mesh,  convergence  of 
the  stress  state  requires  a  more  heavily  refined  mesh  because  of  the  singular  nature 
of  the  stresses  near  the  crack  tip.  The  models  using  the  coarse  mesh  have  approx¬ 
imately  100,000  elements  with  an  element  size  around  the  crack  tip  of  0.025  and 
a  maximum  size  of  0.125.  The  refined  mesh  model  used  for  the  stress  compari¬ 
son  has  approximately  3.9  million  elements  with  the  element  size  varying  between 
0.015  and  0.05. 
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Fig.  15  Schematic  of  the  Abaqus  FE  model 


Abaqus  computes  the  SIFs  for  a  requested  number  of  contours.  The  SIFs  are  typi¬ 
cally  inaccurate  for  the  contours  nearest  to  the  crack  tip  and  often  oscillates  around 
the  correct  value.35,36  For  the  purposes  of  this  work,  the  SIFs  are  calculated  for  10 
contours,  with  the  values  from  the  2  contours  nearest  to  the  crack  tip  being  omitted 
due  to  their  inaccuracy,  and  the  remaining  8  contours  being  averaged  for  the  re¬ 
ported  SIF  value.  For  the  coarse  mesh,  the  SIF  values  for  a/h  =  1.0, 1.25, 1.66  are 
shown  in  Table  2.  The  refined  mesh  is  used  to  calculate  the  SIFs  for  a/ h  =  1.66  us¬ 
ing  the  same  averaging  procedure,  and  Table  3  shows  the  average  values  compared 
with  both  the  coarse  mesh  and  the  values  determined  using  Eq.  28.  Since  the  coarse 
mesh  produced  comparable  results  to  the  solution  using  Eq.  28,  the  normalized  SIFs 
for  a/h  =  1.0, 1.25  are  not  calculated  for  the  refined  mesh.  The  mesh  refinement 
does  have  an  effect  on  the  stress  field  near  the  crack.  The  coarse  mesh  does  not 
produce  smooth  stress  results  for  comparison  to  the  numerical  Gauss-Chebyshev 
solution,  so  the  stress  solution  for  the  refined  mesh  is  used. 
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Table  2  Normalized  SIFs  from  Abaqus  contour  integration 


a/h 


1.0  1.25  1.66 


Contour 

F, 

Fn 

Fi 

Fn 

Fi 

Fn 

1 

1.0454 

-0.0208 

1.0590 

-0.1294 

1.1397 

-0.2048 

2 

0.9535 

-0.1047 

0.9253 

-0.1338 

0.9247 

-0.0816 

3 

0.8689 

-0.0545 

0.8316 

-0.0919 

0.7887 

-0.1149 

4 

0.8801 

-0.0742 

0.8660 

-0.0737 

0.8180 

-0.1023 

5 

0.8971 

-0.0506 

0.8643 

-0.0737 

0.8141 

-0.1036 

6 

0.8886 

-0.0666 

0.8604 

-0.0786 

0.8006 

-0.1106 

7 

0.8672 

-0.0588 

0.8361 

-0.0850 

0.8062 

-0.1082 

8 

0.8886 

-0.0566 

0.8553 

-0.0845 

0.7983 

-0.1126 

9 

0.8790 

-0.0609 

0.8497 

-0.0761 

0.7961 

-0.1126 

10 

0.8762 

-0.0576 

0.8440 

-0.0824 

0.7955 

-0.0990 

Avg  (3-10) 

0.8807 

-0.0600 

0.8510 

-0.0798 

0.8022 

-0.1078 

Table  3  Comparison  of  the  effect  of  Abaqus  mesh  refinement  for  a/h  =  1.66.  The  average 
values  for  Fi  and  f)  /  compared  with  the  Eq.  28  values. 


Coarse  Mesh 

Refined  Mesh 

Numerical  values  Eq.  28a 

F, 

Fn 

F, 

Fn 

Fi  Fn 

0.8022 

-0.1078 

0.8100 

-0.1098 

0.7860  -0.1011 

a  Numerical  values  are  from  Table  B-l,  N  =  9. 
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9.  Conclusions 


In  this  report,  the  method  of  integral  transforms21  and  methods  developed  by  Er- 
dogan  and  Gupta22,23  were  used  to  solve  the  problem  of  a  crack  parallel  to  a  rigid 
boundary  under  remote  tension.  We  derived  a  system  of  singular  integral  equations 
of  the  first  kind,  specific  to  the  problem  at  hand,  which  we  numerically  solved 
using  Gauss-Chebyshev  integration.  SIFs  were  calculated  for  this  problem  which 
were  in  excellent  agreement  with  those  derived  by  others  using  different  solution 
methods.25^27 

Expressions  for  the  normal  stress  ayy(x,y)  and  shear  stress  axy(x,y )  fields  were 
derived  for  the  region  between  the  crack  and  the  rigid  boundary  (shear-stress  free 
symmetry  plane).  We  demonstrated  that  this  region  is  indeed  in  a  state  of  compres¬ 
sion,  as  our  prior  FE  simulations  indicated1  and  despite  the  remotely  applied  tensile 
prestress.  Full-field  elastostatic  stress  solutions  to  stationary  crack  problems  involv¬ 
ing  boundaries,  such  as  those  derived  in  this  report,  are  somewhat  rare  to  find  in  the 
literature,  although  such  solutions  are  necessary  for  verification  of  FE  code  simula¬ 
tions  involving  cracks;  most  papers  of  this  nature  involve  calculation  and  tabulation 
of  only  the  SIFs. 

We  specialized  our  results  in  Eq.  31  and  Eq.  32  to  the  problem  of  a  crack  in  an 
infinite  plate  under  remote  tension  and  showed  that  the  relative  error  in  our  numer¬ 
ically  derived  solutions  are  within  machine  precision  of  the  closed-form  analytical 
solutions  found  in  Zehnder.28  We  also  demonstrated  that  both  the  SIFs  and  stress 
fields  derived  via  numerical  solution  of  the  singular  integral  equations,  compared 
well  with  those  determined  using  the  commercially  available  Abaqus2  FE  code. 
Using  these  results,  we  estimated  the  relative  error  between  the  normal  and  shear 
stress  distributions  using  our  numerical  solution  of  the  singular  integral  equations 
and  Abaqus  FE  results;  it  was  found  that  with  a  relatively  refined  FE  mesh,  Abaqus 
provided  accurate  estimates  of  the  normal  cryy(x,  y)  and  shear  stress  axy(x,  y)  fields 
in  the  vicinity  of  a  stationary  crack  modeled  using  the  XFEM.  Still  requiring  verifi¬ 
cation  with  the  XFEM  in  Abaqus  are  problems  involving  elastodynamic  fields  that 
impinge  upon  stationary  or  dynamically  propagating  cracks. 
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Appendix  A.  Coefficient  Definitions  for  Equation  27 


Approved  for  public  release;  distribution  is  unlimited. 


33 


In  this  appendix,  we  provide  the  additional  equations  needed  to  numerically  evalu¬ 
ate  the  constants  Ak  and  Bp,-  from  Eq.  27.  On  substituting  these  constants  into  the 
definite  integrals  in  Eq.  10  and  Eq.  12,  the  stresses  and  displacements  anywhere 
in  the  domain  can  be  determined  (cf.  Section  5  where  we  calculate  the  normal 
stresses  (Jyy{x,  y),  and  shear  stresses  axy(x,y)  in  the  finite  strip  between  the  crack 
line  y  —  0  and  rigid  boundary  y  =  h  and  compare  our  results  with  computational 
finite  element  results  using  Abaqus). 


The  definitions  of  constants  akn,  bkn,  ckn,  and  dkn  that  appear  in  Eq.  27  can  be  found 
on  page  394  of  Erdogan  et  al.1  are  given  as 


^ kn 


U2k-i{x)Hl1(x)Vl  -  x2  dx  , 


hn  =  /  U2k-i{x)Hl2(x)Vl  -x2dx, 


’-i 


Ckn 


U2k-2(x)H^l1(x)Vl  -  X2  dx  , 


'-i 


(A-l) 


dkn 


U2k-2(x)H^2(x)Vl  -  X2  dx  , 


'-1 


where  U2k-i(x)  and  U2k-2(x)  are  Chebyshev  polynomials  of  the  second  kind,  and 

'Erdogan  F,  Gupta  GD,  Cook  TS.  Numerical  solution  of  singular  integral  equations.  In:  Sih  GC, 
editor.  Mechanics  of  Fracture  I:  Methods  of  Analysis  and  Solutions  of  Crack  Problems;  Netherlands 
(Feyden):  Noordhoff  International  Publishing;  1973.  p.  368-425. 
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polynomial  functions  Hlj{x)  are  defined  on  page  393  of  Erdogan  et  al.1  as 


-I 


Hn(x)  =  -  j  ku(x,t)T2n(t)  (l  -  t2)  2  dt, 


H„2(x)  =  ~  J  k12(x,t)T2n-i(t)  (l  -  t2)  2  dt, 


(A-2) 


If1  _i 

H^{x)  =  -  J  k2i(x,  t)T2n(t)  (l  —  t2)  2  dt, 


If'1  _i 

=  -  j  k22(x,t)T2n-i(t)  (l  -  t2)  2  dt, 


where  T2n(x )  and  T2n-\  (x)  are  Chebyshev  polynomials  of  the  first  kind,  and  the 
kernel  functions  kij(x,t )  in  Eq.  A-2  are  given  in  Eq.  25.  In  practice,  polynomial 
functions  H^(x)  in  Eq.  A-2  are  determined  in  Mathematica2  using  Gauss-Chebyshev 
numerical  integration  with  n  =  300  terms,3  viz., 

4==  dt  «  ^  vjJ(x,  U)  ,  (A- 3) 

^  *  i=i 

with  weights,  Wi  =  ^  and  abscissas,  f*  =  cos(  )  obtained  from  the  zeros  of 
the  orthogonal  Chebyshev  polynomials  of  the  first  kind  (cf.  Eq.  25.4.38  on  page 
889  of  Abramovitz  and  Stegun4). 

The  polynomial  functions  H'^(x)  determined  with  Eq.  A-3  are  then  substituted 
into  Eq.  A-l  for  determining  constants  a/,„,  ckn,  and  dt,,  used  in  system  Eq.  27. 
Gauss-Chebyshev  numerical  integration  with  n  =  300  terms,  is  again  employed 

2  Mathematica  edition  Ver.  10.0.  Champaign  (IL):  Wolfram  Research;  2013. 

^Numerical  integration  using  n  =  300  terms  brings  the  numerical  solution  accuracy  to  15  sig¬ 
nificant  digits  which  is  at  machine  precision,  cf.  Fig.  10. 

4Abramowitz  M,  Stegun  I.  Handbook  of  mathematical  functions.  9th  ed.  New  York  (NY): 
Dover  Publications;  1970. 
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using  Mathematica2  to  solve  Eq.  A-l  with  approximation, 

/I  n 

f(x,t)V  1  -  t2dt  PH  y ^Wif(x,tj)  ,  (A-4) 

1  i= i 

with  weights,  Wi  =  ^j-j-  sin2(^j-)  and  abscissas,  U  =  cos(^)  obtained  from  the 
zeros  of  the  orthogonal  Chebyshev  polynomials  of  the  second  kind  (cf.  Eq.  25.4.40 
on  page  889  of  Abramovitz  and  Stegun  (1970)4). 

As  an  example,  for  N  =  3  system  Eq.  27  appears  in  matrix-vector  form  as, 


ax  =  b 

where, 


+  All 

A 12 

fl'13 

^11 

012 

O13 

®21 

§  +  ®22 

®23 

&21 

b-22 

b>23 

A31 

®32 

f  +  a33 

^31 

b'32 

b'33 

Cll 

Cl2 

Cl3 

2  ^ 

d\2 

d\3 

C21 

C22 

C23 

^21 

§  +  d-22 

d'23 

C31 

C32 

C33 

^31 

d%2 

?  +  d_ 

"  A1  ' 

0 

A-2 

0 

A.% 

0 

,  b  = 

7T  p  K+l 

2  B 1  2At 

Bi 

7 

b2 

0 

b3 

0 

(A- 5) 


(A-6) 


(A- 7) 
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The  vector  x  or  the  constants  A{.  A2,  A3,  B\ .  B2,  B3  for  iV  =  3  are  obtained  by 
finding  the  inverse  of  non-symmetric  constant  matrix  a  given  by  Eq.  A-6,  and  post- 
multiplying  this  by  b  defined  by  F\ and  F-2k  in  Eq.  27,  i.e.,  x  =  a-1  b.  Given  the 
Ai,  A2,  A3,  B\ ,  B2,  B3  enables  determination  of  approximate  values  for  the  disloca¬ 
tion  densities  in  Eq.  26,  which  together  with  A,3  (£)  in  Eq.  18  enables  determination 
of  the  displacement  and  stresses  throughout  the  medium.  In  addition,  the  calculation 
of  stress  intensity  factors  (SIFs)  is  rather  straightforward  and  discussed  in  Section  4. 
Table  B-l  in  Appendix  B  illustrates  SIF  convergence  for  N  =  3,  6,  9. 
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Appendix  B.  Stress  Intensity  Factor  Convergence 
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Table  B-l  Normalized  SIFs  for  2  parallel  cracks  in  an  infinite  plane  under  remote  tension, 
where  Kj  =  Fj a  \/Fo,  Kn  =  Fua^pna,  and  k\  =  Ff;  k2  =  Fn  are  determined  using  Eq.  28 


a/h 

Fj  ( N  =  3) 

Fi  {N  =  6) 

F/  (N  =  9) 

CO 

II 

co' 

II 

Fn  (N  =  9) 

0.01 

0.99996 

0.99996 

0.99996 

-1.8747xl0~7 

-1.8747  xl0~7 

-1.8747x  10~7 

0.1 

0.9963 

0.9963 

0.9963 

-0.00018 

-0.00018 

-0.00018 

0.2 

0.9857 

0.9857 

0.9857 

-0.0014 

-0.0014 

-0.0014 

0.4 

0.9505 

0.9505 

0.9505 

-0.0094 

-0.0094 

-0.0094 

0.6 

0.9086 

0.9086 

0.9086 

-0.0246 

-0.0246 

-0.0246 

0.8 

0.8706 

0.8707 

0.8707 

-0.0429 

-0.0429 

-0.0429 

1.0 

0.8403 

0.8404 

0.8404 

-0.0607 

-0.0607 

-0.0607 

1.25 

0.8130 

0.8131 

0.8131 

-0.0793 

-0.0793 

-0.0793 

1.66 

0.7864 

0.7860 

0.7860 

-0.1008 

-0.1011 

-0.1011 

2.0 

0.7754 

0.7737 

0.7737 

-0.1122 

-0.1128 

-0.1128 

5.0 

0.7504 

0.7441 

0.7445 

-0.1649 

-0.1566 

-0.1567 

10.0 

0.6941 

0.7441 

0.7385 

-0.1992 

-0.1785 

-0.1792 
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